Educational notebook based on Sargolini 2006¶
This tutorial demonstrates how to access the dataset published in Sargolini et al. (Science, 2006) using dandi.
The dataset contains spike times for recorded grid cells from the medial entorhinal cortex (MEC) in rats that explored two-dimensional environments. The behavioral data includes position from the tracking LED(s).
Contents:¶
Streaming NWB files ¶
This section demonstrates how to access the files on DANDI without downloading them.
Based on the Streaming NWB files tutorial from PyNWB.
The DandiAPIClient can be used to get the S3 URL of the NWB file stored in the DANDI Archive.
from dandi.dandiapi import DandiAPIClient
dandiset_id, nwbfile_path = "000582", "sub-10073/sub-10073_ses-17010302_behavior+ecephys.nwb" # file size ~15.6MB
# Get the location of the file on DANDI
with DandiAPIClient() as client:
asset = client.get_dandiset(dandiset_id, 'draft').get_asset_by_path(nwbfile_path)
s3_url = asset.get_content_url(follow_redirects=1, strip_query=True)
Create a virtual filesystem using fsspec which will take care of requesting data from the S3 bucket whenever data is read from the virtual file.
from fsspec.implementations.cached import CachingFileSystem
from fsspec import filesystem
from h5py import File
from pynwb import NWBHDF5IO
# first, create a virtual filesystem based on the http protocol
fs=filesystem("http")
# create a cache to save downloaded data to disk (optional)
fs = CachingFileSystem(
fs=fs,
cache_storage="nwb-cache", # Local folder for the cache
)
file_system = fs.open(s3_url, "rb")
file = File(file_system, mode="r")
# Open the file with NWBHDF5IO
io = NWBHDF5IO(file=file, load_namespaces=True)
nwbfile = io.read()
nwbfile
root (NWBFile)
file_create_date
acquisition (1)
ElectricalSeries
data
electrodes
table
id
keywords
processing (2)
behavior
data_interfaces (1)
Position
spatial_series (1)
SpatialSeriesLED1
data
timestamps
ecephys
data_interfaces (1)
LFP
electrical_series (1)
ElectricalSeriesLFP
data
electrodes
table
id
electrodes
id
electrode_groups (1)
ElectrodeGroup
device
devices (1)
EEG
subject
units
id
Access data and metadata ¶
This section demonstrates how to access the data in the NWB files.
Subject¶
The nwbfile.subject field holds information about the experimental subject, such as age (in ISO 8601 Duration format), sex, and species in latin binomial nomenclature.
nwbfile.subject
subject (Subject)
Accessing behavior data ¶
The "behavior" processing module holds the behavior data in the NWB file which can be accessed as
nwbfile.processing["behavior"].
Position¶
The position data is stored in a SpatialSeries object which can be accessed from the "Position" container as nwbfile.processing["behavior"]["Position"].
Note that not all sessions have position data from two tracking LEDs.
spatial_series = nwbfile.processing["behavior"]["Position"]["SpatialSeriesLED1"]
spatial_series
SpatialSeriesLED1 (SpatialSeries)
data
timestamps
# Showing the trace of position
import plotly.express as px
fig = px.line(
x=spatial_series.data[:, 0] * spatial_series.conversion,
y=spatial_series.data[:, 1] * spatial_series.conversion,
labels=dict(x="X (meters)", y="Y (meters)"),
)
fig.update_layout(yaxis=dict(scaleanchor="x"))
fig.show()
Accessing spike times ¶
The Units table holds the spike times which can be accessed as nwbfile.units and can also be converted to a pandas DataFrame.
nwbfile.units.to_dataframe()
| unit_name | spike_times | histology | hemisphere | depth | |
|---|---|---|---|---|---|
| id | |||||
| 0 | t1c1 | [0.7903958333333333, 0.794, 0.8111666666666667... | MEC LII | 0.0024 | |
| 1 | t2c1 | [1.0451354166666667, 1.7003854166666668, 2.315... | MEC LII | 0.0024 | |
| 2 | t2c3 | [0.18273958333333334, 0.5340729166666667, 0.57... | MEC LII | 0.0024 | |
| 3 | t3c1 | [1.0358229166666666, 1.04803125, 1.69642708333... | MEC LII | 0.0024 | |
| 4 | t3c2 | [2.43025, 2.4398333333333335, 3.17965625, 3.39... | MEC LII | 0.0024 | |
| 5 | t3c3 | [2.1157708333333334, 2.425427083333333, 3.3630... | MEC LII | 0.0024 | |
| 6 | t3c4 | [0.07945833333333334, 2.244947916666667, 3.173... | MEC LII | 0.0024 | |
| 7 | t4c1 | [2.4301666666666666, 2.439770833333333, 3.1795... | MEC LII | 0.0024 |
Visualizing rate maps ¶
This section demonstrates how to show the rate maps of the recorded cells.
We will use PYthon Neural Analysis Package (pynapple) to calculate the rate maps.
Using the compute_2d_tuning_curves() function we can compute firing rate as a function of position (map of neural activity as the animal explored the environment).
#!pip install pynapple
import pynapple as nap
position_over_time = nap.TsdFrame(
d=spatial_series.data[:],
t=spatial_series.timestamps[:],
columns=["x","y"],
)
spike_times_group = nap.TsGroup({cell_id: nap.Ts(spikes) for cell_id, spikes in enumerate(nwbfile.units["spike_times"])})
num_bins = 15
rate_maps, position_bins = nap.compute_2d_tuning_curves(
spike_times_group,
position_over_time,
num_bins,
)
import plotly.express as px
cell_ind = 1
unit_name = nwbfile.units["unit_name"][cell_ind]
layer = nwbfile.units["histology"][cell_ind]
fig = px.imshow(rate_maps[cell_ind], color_continuous_scale="viridis")
fig.update_xaxes(showticklabels=False)
fig.update_yaxes(showticklabels=False)
fig.update_layout(title=f"{unit_name} from {layer}", coloraxis_showscale=False)
fig.show()

Visualizing grid cells activity¶
To determine whether the firing fields of individual cells formed a grid structure, we will calculate the spatial autocorrelation for the rate map of each cell.
The autocorrelograms are based on Pearson’s product moment correlation coefficient with corrections for edge effects and unvisited locations. With λ (x, y) denoting the average rate of a cell at location (x, y), the autocorrelation between the fields with spatial lags of τx and τy was estimated as:

where the summation is over all n pixels in λ (x, y) for which rate was estimated for both λ (x, y) and λ (x - τx, y - τy). Autocorrelations were not estimated for lags of τx, τy where n < 20.
The degree of spatial periodicity (gridness) can be determined for each cell by rotating the autocorrelation map for each cell in steps of 6 degrees (from 0 to 180 degrees) and computing the correlation between the rotated map and the original. The correlation is confined to the area defined by a circle around the peaks that are closest to the centre of the map, and the central peak is not included in the analysis.
The ‘gridness’ of a cell can be expressed as the difference between the lowest correlation at 60 and 120 degrees (where a peak correlation would be expected due to the triangular nature of the grid) and the highest correlation at 30, 90, and 150 degrees (where the minimum correlation would be expected). When the correlations at 60 and 120 degrees of rotation exceeded each of the correlations at 30, 90 and 150 degrees (gridness > 0), the cell was classified as a grid cell.

import numpy as np
def create_coer_arr(arr, rad_min=None, rad_max=None):
""" Creating an array for correlation(tau_x, tau_y)
I will take tau_x from the range (-arr.shape[0]+1, arr.shape[0]-1) and the same for tau_y
"""
sh_x, sh_y = arr.shape
# creating an array full of nan's
coer_arr = np.full((2*sh_x-1, 2*sh_y-1), np.nan)
for ii in range(0, 2*(sh_x-1)):
for jj in range(0, 2*(sh_y-1)):
# shifting tau_x/y
tau_x = ii-sh_x+1
tau_y = jj-sh_y+1
# if rad_max and rad_min is provided, I only calculate the correlation for points between rad_min and rad_max
if rad_max is not None and ((tau_x**2 + tau_y**2)**0.5 > rad_max):
continue
if rad_min is not None and ((tau_x**2 + tau_y**2)**0.5 < rad_min):
continue
coer_arr[ii, jj] = pearson_autocor(arr, lag_x=tau_x, lag_y=tau_y)
return coer_arr
def pearson_autocor(arr, lag_x, lag_y):
""" Calculating Pearson autocorrelation for an array that can contain NaN values."""
sh_x, sh_y = arr.shape
if abs(lag_x) >= sh_x or abs(lag_y) >= sh_y:
raise Exception(f"abs(lag_x), abs(lag_y) have to be smaller than {sh_x}, {sh_y}, but {lag_x}, {lag_y} provided")
# calculating sum for elements that meet the requirements
n = 0
sum1, sum2, sum3, sum4, sum5 = 0, 0, 0, 0, 0
for ii in range(0, sh_x):
for jj in range(0, sh_y):
# checking if the indices are within the array
if 0 <= ii-lag_x < sh_x and 0 <= jj-lag_y < sh_y:
# checking if both values (in ii,jj and shifted) are not nan
if not np.isnan(arr[ii, jj]) and not np.isnan(arr[ii-lag_x, jj-lag_y]):
n += 1
sum1 += arr[ii, jj] * arr[ii-lag_x, jj-lag_y]
sum2 += arr[ii, jj]
sum3 += arr[ii-lag_x, jj-lag_y]
sum4 += (arr[ii, jj])**2
sum5 += (arr[ii-lag_x, jj-lag_y])**2
# according to the paper they had this limit for number of points
if n < 20:
return np.nan
numerator = n * sum1 - sum2 * sum3
denominator = (n * sum4 - sum2**2)**0.5 * (n * sum5 - sum3**2)**0.5
cor = numerator / denominator
return cor
import plotly.express as px
from plotly.subplots import make_subplots
import plotly.graph_objects as go
from scipy.ndimage import rotate
import scipy
rate_maps_50_bin, _ = nap.compute_2d_tuning_curves(
spike_times_group,
position_over_time,
50,
)
for cell_ind in range(len(rate_maps)):
unit_name = nwbfile.units["unit_name"][cell_ind]
layer = nwbfile.units["histology"][cell_ind]
fig = make_subplots(
rows=1,
cols=3,
subplot_titles=(f'{unit_name} rate map', f'{unit_name} auto-correlation', "periodicity"),
)
rate_map_plot = go.Heatmap(z=rate_maps[cell_ind], colorscale='Viridis', showscale=False)
fig.add_trace(rate_map_plot, row=1, col=1)
# Compute auto-correlation
autocorr = create_coer_arr(rate_maps_50_bin[cell_ind], rad_max=34, rad_min=6)
autocorr_nonan = np.nan_to_num(autocorr, copy=True, nan=0.0)
correlations = []
angles = np.arange(0, 186, 6)
for angle in angles:
autocorr_nonan_rotated = rotate(autocorr_nonan, angle=angle, reshape=False)
cor = scipy.stats.pearsonr(autocorr_nonan_rotated.flatten(), autocorr_nonan.flatten())[0]
correlations.append(cor)
gridness = max(correlations[10], correlations[20]) - max(correlations[5], correlations[15], correlations[25])
gridness = np.round(gridness, 2)
autocorr_rate_map = go.Heatmap(z=autocorr, colorscale='Viridis', showscale=False)
fig.add_trace(autocorr_rate_map, row=1, col=2)
line_trace = go.Scatter(
x=angles,
y=correlations,
mode='lines',
marker=dict(color="black"),
)
fig.add_trace(line_trace, row=1, col=3)
fig.update_xaxes(showticklabels=False)
fig.update_yaxes(showticklabels=False)
fig.update_xaxes(showticklabels=True, row=1, col=3, title_text="Rotation (deg)")
fig.update_yaxes(showticklabels=True, row=1, col=3, title_text="Correlation (r)")
fig.update_layout(
title=f"{layer}: Gridness score {gridness}",
xaxis3 = dict(
tickmode="array",
tickvals=[30, 60, 90, 120, 150, 180],
)
)
fig.show()
/Users/weian/anaconda3/envs/turner311/lib/python3.11/site-packages/pynapple/process/tuning_curves.py:202: RuntimeWarning: invalid value encountered in divide

